/*
Copyright (c) 2009-2010, Warwick Warp Limited
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of the Warwick Warp Limited nor the
      names of its contributors may be used to endorse or promote products
      derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL WARWICK WARP LIMITED BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#ifndef WARRAY_WARRAY_BLAS_H_INCLUDED
#define WARRAY_WARRAY_BLAS_H_INCLUDED

#include <complex>
#include <numeric>
#include "warray_help.h"

namespace warp {
namespace blas {
namespace detail {
template<typename T>
struct auto_return_type { typedef T type; };
template<typename T>
struct auto_return_type< std::complex<T> > 
{ 
	typedef std::complex<typename auto_return_type<T>::type> type;
};
template<> struct auto_return_type<int> { typedef double type; };
template<> struct auto_return_type<unsigned int> { typedef double type; }; 
template<> struct auto_return_type<short> { typedef double type; };
template<> struct auto_return_type<unsigned short> { typedef double type; };
template<> struct auto_return_type<signed char> { typedef double type; };
template<> struct auto_return_type<unsigned char> { typedef double type; };
} // namespace detail

enum conjtrans_e
{
	dont_transpose=111,
	transpose=112,
	conjugate_transpose=113
};

enum uplo_e
{
	upper_triangular=121,
	lower_triangular=122
};

// This argument type indicates that conjugate_transpose isn't accepted
typedef conjtrans_e trans_e;

// Level 1 routines
//! \brief Compute \f$ \sum_{i=0}^{n-1} \left|x_i\right| \f$
template<typename T>
T asum(size_t n, const T *x, int incx)
{
	T result = T();
	for(;n>0;--n,x+=incx) 
		result += std::abs(*x);
	return result;
}

//! \brief Compute \f$ \sum_{i=0}^{n-1} \left|\text{real}(x_i)\right| + \left|\text{imag}(x_i)\right| \f$
template<typename T>
T asum(size_t n, const std::complex<T> *x, int incx)
{
	T result = T();
	for(;n>0;--n,x+=incx) 
		result += std::abs(x->real()) + std::abs(x->imag());
	return result;
}

//! \brief Compute \f$ y_i = y_i + ax_i \f$
template<typename T>
void axpy(size_t n, const T *x, int incx, T *y, int incy, T a=T(1))
{
	for(;n>0;--n,x+=incx,y+=incy)
		*y += a**x;
}

//! \brief Compute \f$ y_i = x_i \f$
template<typename T>
void copy(size_t n, const T *x, int incx, T *y, int incy)
{
	for(;n>0;--n,x+=incx,y+=incy)
		*y = *x;
}

//! \brief Compute \f$ \sum_{i=0}^{n-1} x_iy_i \f$
template<typename T>
T dot(size_t n, const T *x, int incx, const T *y, int incy)
{
	T result = T();
	for(;n>0;--n,x+=incx,y+=incy)
		result += *x * *y;
	return result;
}

//! \brief Compute \f$ \sum_{i=0}^{n-1} \text{conj}(x_i)y_i \f$
template<typename T>
std::complex<T> dotc(size_t n, const std::complex<T> *x, int incx, const std::complex<T> *y, int incy)
{
	std::complex<T> result = std::complex<T>();
	for(;n>0;--n,x+=incx,y+=incy)
		result += std::conj(*x) * *y;
	return result;
}

//! \brief Compute \f$ \sqrt{\sum_{i=0}^{n-1} x_i^2} \f$
template<typename T>
typename detail::auto_return_type<T>::type nrm2(size_t n, const T *x, int incx)
{
	typename detail::auto_return_type<T>::type dp( dot(n,x,incx,x,incx) );
	return std::sqrt(dp);
}

//! \brief Compute \f$ \sqrt{\sum_{i=0}^{n-1} \text{real}(x_i)^2 + \text{imag}(x_i)^2} \f$
template<typename T>
typename detail::auto_return_type<T>::type nrm2(size_t n, const std::complex<T> *x, int incx)
{
	typename detail::auto_return_type<T>::type result = T();
	for(;n>0;--n,x+=incx)
		result += x->real()*x->real() + x->imag()*x->imag();
	return std::sqrt(result);
}

/*!
	\brief Compute
	\f$ x_i^\prime = c x_i + s y_i \f$
	\f$ y_i^\prime = c y_i - s x_i \f$
	
	Note that this is a CLOCKWISE rotation of the vector (x,y), which is the 
	opposite of our standard anticlockwise rotation specifications.
*/
template<typename T>
void rot(size_t n, T *x, int incx, T *y, int incy, const T& c, const T& s)
{
	for(;n>0;--n,x+=incx,y+=incy)
	{
		T xx = c**x + s**y;
		T yy = c**y - s**x;
		*x = xx;
		*y = yy;
	}
}

//! \brief Compute \f$ x_i = ax_i \f$
template<typename T>
void scal(size_t n, const T& a, T *x, int incx)
{
	for(;n>0;--n,x+=incx)
		*x *= a;
}

//! \brief Compute \f$ x_i = ax_i \f$
template<typename T>
void scal(size_t n, const T& a, std::complex<T> *x, int incx)
{
	for(;n>0;--n,x+=incx)
		*x *= a;
}


/*!
	\brief Compute
	\f$ x_i^\prime = y_i \f$
	\f$ y_i^\prime = x_i \f$
*/
template<typename T>
void swap(size_t n, T *x, int incx, T *y, int incy)
{
	for(;n>0;--n,x+=incx,y+=incy)
		std::swap(*x,*y);
}

//! \brief Compute \text{argmax}_i \left|x_i\right|
template<typename T>
ptrdiff_t amax(size_t n, const T *x, int incx)
{
	if(n<=1)
		return 0;
	ptrdiff_t result = 0;
	T         xmax   = std::abs(*x);
	ptrdiff_t index  = 0;
	for(--n,++index,x+=incx;n>0;--n,++index,x+=incx)
	{
		T xabs = std::abs(*x);
		if(xabs > xmax)
		{
			xmax   = xabs;
			result = index;
		}
	}
	return result;
}

//! \brief Compute \text{argmax}_i \left|\text{real}(x_i)\right| + \left|\text{imag}(x_i)\right|
template<typename T>
ptrdiff_t amax(size_t n, const std::complex<T> *x, int incx)
{
	if(n<=1)
		return 0;
	ptrdiff_t result = 0;
	T         xmax   = std::abs(x->real())+std::abs(x->imag());
	ptrdiff_t index  = 0;
	for(--n,++index,x+=incx;n>0;--n,++index,x+=incx)
	{
		T xabs = std::abs(x->real())+std::abs(x->imag());
		if(xabs > xmax)
		{
			xmax   = xabs;
			result = index;
		}
	}
	return result;
}

//! \brief Compute \text{argmin}_i \left|x_i\right|
template<typename T>
ptrdiff_t amin(size_t n, const T *x, int incx)
{
	if(n<=1)
		return 0;
	ptrdiff_t result = 0;
	T         xmin   = std::abs(*x);
	ptrdiff_t index  = 0;
	for(--n,++index,x+=incx;n>0;--n,++index,x+=incx)
	{
		T xabs = std::abs(*x);
		if(xabs < xmin)
		{
			xmin   = xabs;
			result = index;
		}
	}
	return result;
}

//! \brief Compute \text{argmin}_i \left|\text{real}(x_i)\right| + \left|\text{imag}(x_i)\right|
template<typename T>
ptrdiff_t amin(size_t n, const std::complex<T> *x, int incx)
{
	if(n<=1)
		return 0;
	ptrdiff_t result = 0;
	T         xmin   = std::abs(x->real())+std::abs(x->imag());
	ptrdiff_t index  = 0;
	for(--n,++index,x+=incx;n>0;--n,++index,x+=incx)
	{
		T xabs = std::abs(x->real())+std::abs(x->imag());
		if(xabs < xmin)
		{
			xmin   = xabs;
			result = index;
		}
	}
	return result;
}

// Level 2 routines
/*!
	\brief General matrix vector product.
	
	Compute \f$ y = \alphaAx + \betay \f$ or \f$ y = \alpha^\text{T}Ax + \betay \f$
	depending on the value of the transpose flag.
	
	The dimensions m and n specify the dimensions of the matrix A regardless of 
	the value of the transpose parameter.
	
	This implementation does NOT support aliasing of x and y or a and y.
*/
template<typename T>
void gemv( size_t m, size_t n, const T* a, int lda, const T* x, int incx, T* y, int incy, const T& alpha=T(1), const T& beta=T(0), trans_e trans=dont_transpose)
{
	int inca = trans!=dont_transpose ? lda : 1;
	int rowa = trans!=dont_transpose ? 1 : lda;
	size_t rows = trans!=dont_transpose ? n : m;
	size_t cols = trans!=dont_transpose ? m : n;
	warp::detail::alias_check(x,incx,cols,y,incy,rows);
	for(;rows>0;--rows,a+=rowa,y+=incy)
		*y = alpha*dot(cols,a,inca,x,incx) + beta**y;
}

/*!
	\brief General rank 1 matrix update.
	
	Compute \f$ A = A + \alphaxy^\text{T}
	
	This implementation does not support aliasing of a and x or y.
*/
template<typename T>
void ger(size_t m, size_t n, const T* x, int incx, const T* y, int incy, T* a, int lda, const T& alpha=T(1))
{
	const T* xi=x;
	T* ai=a;
	for(size_t i=0;i<m;++i,xi+=incx,ai+=lda)
	{
		const T* yj=y;
		T* aij=ai;
		for(size_t j=0;j<n;++j,yj+=incy,++aij)
		{
			*aij += alpha * *xi * *yj;
		}
	}
}

/*!
	\brief General conjugated rank 1 matrix update.
	
	Compute \f$ A = A + \alphax\text{conj}(y)^\text{T}
	
	This implementation does not support aliasing of a and x or y.
*/

template<typename T>
void gerc(size_t m, size_t n, const std::complex<T>* x, int incx, const std::complex<T>* y, int incy, std::complex<T>* a, int lda,  const std::complex<T>& alpha=std::complex<T>(T(1)))
{
	const std::complex<T>* xi=x;
	std::complex<T>* ai=a;
	for(size_t i=0;i<m;++i,xi+=incx,ai+=lda)
	{
		const std::complex<T>* yj=y;
		std::complex<T>* aij=ai;
		for(size_t j=0;j<n;++j,yj+=incy,++aij)
		{
			*aij += alpha * *xi * std::conj(*yj);
		}
	}
}

/*!
	\brief Symmetric matrix vector product.
	
	Compute \f$ y = \alphaAx + \betay \f$.
	
	This implementation does NOT support aliasing of x and y or a and y.
	
	Note that this implementation is slower than the general matrix vector 
	product because it must perform the same calculation, but track two pointers 
	into the array a to ensure that only the upper (or lower) triangular part is 
	accessed.
*/

template<typename T>
void symv(size_t n, const T* a, int lda, const T* x, int incx, T* y, int incy, const T& alpha=T(1), const T& beta=T(0), uplo_e upper=upper_triangular)
{
	warp::detail::alias_check(x,incx,n,y,incy,n);
	
	const T *a0i=a; // Pointer to the start of the ith column
	const T *ai0=a; // Pointer to the start of the ith row
		
	for(size_t i=0;i<n;++i,++a0i,ai0+=lda,y+=incy)
	{
		T val = T(0);
		const T* aij=ai0; // Pointer to the element j,i
		const T* aji=a0i; // Pointer to the element i,j
		const T* xi=x;
		
		// Lower triangle (j<i)
		for(size_t j=0;j<i;++j,aji+=lda,++aij,xi+=incx)
			val += (upper==upper_triangular ? *aji : *aij) * *xi;
		// Upper triangle (j>=i)
		for(size_t j=i;j<n;++j,aji+=lda,++aij,xi+=incx)
			val += (upper==upper_triangular ? *aij : *aji) * *xi;
		
		*y = alpha*val + beta**y;
	}
}

/*!
	\brief Perform a rank 1 update of a symmetric matrix
	
	Compute \f$ A = A + \alphaxx^\text{T}
	
	This implementation does not support aliasing of a and x.
*/
template<typename T>
void syr(size_t n, const T *x, int incx, T *a, int lda, const T& alpha = T(1), uplo_e upper=upper_triangular)
{
	if(upper==upper_triangular)
	{
		const T* xi=x;
		T* ai=a;
		for(size_t i=0;i<n;++i,xi+=incx,ai+=lda)
		{
			const T* xj=xi;
			T* aij=ai+i;
			for(size_t j=i;j<n;++j,xj+=incx,++aij)
			{
				*aij += alpha * *xi * *xj;
			}
		}
	}
	else
	{
		const T* xi=x;
		T* ai=a;
		for(size_t i=0;i<n;++i,xi+=incx,ai+=lda)
		{
			const T* xj=x;
			T* aij=ai;
			for(size_t j=0;j<=i;++j,xj+=incx,++aij)
			{
				*aij += alpha * *xi * *xj;
			}
		}
	}
}

/*!
	\brief Perform a rank 2 update of a symmetric matrix
	
	Compute \f$ A = A + \alphaxy^\text{T} + \alphayx^\text{T}
	
	This implementation does not support aliasing of a and x or a and y.
*/

template<typename T>
void syr2(size_t n, const T *x, int incx, const T *y, int incy, T *a, int lda, const T& alpha=T(1), uplo_e upper=upper_triangular)
{
	if(upper==upper_triangular)
	{
		const T* xi=x;
		const T* yi=y;
		T* ai=a;
		for(size_t i=0;i<n;++i,xi+=incx,yi+=incy,ai+=lda)
		{
			const T* xj=xi;
			const T* yj=yi;
			T* aij=ai+i;
			for(size_t j=i;j<n;++j,xj+=incx,yj+=incy,++aij)
			{
				*aij += alpha * (*xi * *yj + *xj * *yi);
			}
		}
	}
	else
	{
		const T* xi=x;
		const T* yi=y;
		T* ai=a;
		for(size_t i=0;i<n;++i,xi+=incx,yi+=incy,ai+=lda)
		{
			const T* xj=x;
			const T* yj=y;
			T* aij=ai;
			for(size_t j=0;j<=i;++j,xj+=incx,yj+=incy,++aij)
			{
				*aij += alpha * (*xi * *yj + *xj * *yi);
			}
		}
	}
}

// Level 3 routines

/*!
	\brief Perform a general matrix-matrix multiply
	
	\f$ C = \alphaAB + \betaC \f$
	
	A (or A' if transa is true) is an m by k matrix
	B (or B' if transb is true) is a k by n matrix
	C is an m by n matrix
	
	This implementation does not support aliasing of c with either of the input 
	arrays a or b.
	
	NB: conjugate-transpose isn't implemented (yet)
*/
template<typename T>
void gemm(size_t m, size_t n, size_t k, const T *a, int lda, const T *b, int ldb, T *c, int ldc, const T& alpha=T(1), const T& beta=T(0), conjtrans_e transa=dont_transpose, conjtrans_e transb=dont_transpose)
{
	for(;n>0;--n,++c,b+=(transb==transpose?ldb:1))
	{
		// Compute each column of c by computing C_i = AB_i
		gemv(transa==transpose?k:m,transa==transpose?m:k, // Ug, the BLAS routines are inconsistent in their handling of transposition
			 a,lda,
			 b,transb==transpose?1:ldb,
			 c,ldc,
			 alpha,beta,
			 transa);
		
	}
}

/*!
	\brief Perform a rank-k update of a symmetric matrix
	
	\f$ C = \alpha AA^\prime + \beta C \f$
	or
	\f$ C = \alpha A^\primeA + \beta C \f$
	
	A is an n by k matrix if trans is false and a k by n matrix if trans is true.
	
*/
template<typename T>
void syrk(size_t n, size_t k, const T *a, int lda, T *c, int ldc, const T& alpha=T(1), const T& beta=T(0), trans_e trans=dont_transpose, uplo_e upper=upper_triangular)
{
	T*ci=c;
	int cistep = ldc+1;
	int cjstep = (upper==upper_triangular) ? 1 : ldc;
	
	const T* ai = a;
	int astep = (trans==dont_transpose) ? lda : 1;
	int ainc  = (trans==dont_transpose) ? 1 : lda;
	
	for(size_t i=0;i<n;++i,ci+=cistep,ai+=astep)
	{
		T*cij=ci;
		const T* aj=ai;
		for(size_t j=i;j<n;++j,cij+=cjstep,aj+=astep)
		{
			*cij = alpha * dot(k,ai,ainc,aj,ainc) + beta * *cij;
		}
	}
}

/*!
	\brief Perform a rank-2k update of a symmetric matrix

	\f$ C = \alpha A B^\prime + \alpha B A^\prime + \beta C \f$
	or
	\f$ C = \alpha A^\prime B + \alpha B^\prime A + \beta C \f$
	
	A and B are n by k if trans is false and k by n if trans is true.
*/
template<typename T>
void syr2k(size_t n, size_t k, const T *a, int lda,const T *b, int ldb, T *c, int ldc, const T& alpha=T(1), const T& beta=T(0), trans_e trans=dont_transpose, uplo_e upper=upper_triangular)
{
	T*ci=c;
	int cistep = ldc+1;
	int cjstep = (upper==upper_triangular) ? 1 : ldc;
	
	const T* ai = a;
	const T* bi = b;
	int astep = (trans==dont_transpose) ? lda : 1;
	int ainc  = (trans==dont_transpose) ? 1 : lda;
	int bstep = (trans==dont_transpose) ? ldb : 1;
	int binc  = (trans==dont_transpose) ? 1 : ldb;
	
	for(size_t i=0;i<n;++i,ci+=cistep,ai+=astep,bi+=bstep)
	{
		T*cij=ci;
		const T* aj=ai;
		const T* bj=bi;
		for(size_t j=i;j<n;++j,cij+=cjstep,aj+=astep,bj+=bstep)
		{
			*cij = alpha * (dot(k,ai,ainc,bj,binc) + dot(k,bi,binc,aj,ainc)) + beta * *cij;
		}
	}
}

#if defined(USE_ACCELERATE) || defined(USE_MKL) || defined(USE_CBLAS)
template<> float asum<float>(size_t n, const float *x, int incx);
template<> double asum<double>(size_t n, const double *x, int incx);
template<> float asum<float>(size_t n, const std::complex<float> *x, int incx);
template<> double asum<double>(size_t n, const std::complex<double> *x, int incx);
template<> void copy<float>(size_t n, const float *x, int incx, float *y, int incy);
template<> void copy<double>(size_t n, const double *x, int incx, double *y, int incy);
template<> void copy< std::complex<float> >(size_t n, const std::complex<float> *x, int incx, std::complex<float> *y, int incy);
template<> void copy< std::complex<double> >(size_t n, const std::complex<double> *x, int incx, std::complex<double> *y, int incy);
template<> float dot<float>(size_t n, const float *x, int incx, const float *y, int incy);
template<> double dot<double>(size_t n, const double *x, int incx, const double *y, int incy);
template<> std::complex<float> dot< std::complex<float> >(size_t n, const std::complex<float> *x, int incx, const std::complex<float> *y, int incy);
template<> std::complex<double> dot< std::complex<double> >(size_t n, const std::complex<double> *x, int incx, const std::complex<double> *y, int incy);
template<> std::complex<float> dotc<float>(size_t n, const std::complex<float> *x, int incx, const std::complex<float> *y, int incy);
template<> std::complex<double> dotc<double>(size_t n, const std::complex<double> *x, int incx, const std::complex<double> *y, int incy);
template<> detail::auto_return_type<float>::type nrm2<float>(size_t n, const float *x, int incx);
template<> detail::auto_return_type<double>::type nrm2<double>(size_t n, const double *x, int incx);
template<> detail::auto_return_type<float>::type nrm2<float>(size_t n, const std::complex<float> *x, int incx);
template<> detail::auto_return_type<double>::type nrm2<double>(size_t n, const std::complex<double> *x, int incx);
template<> void rot<float>(size_t n, float *x, int incx, float *y, int incy, const float& c, const float& s);
template<> void rot<double>(size_t n, double *x, int incx, double *y, int incy, const double& c, const double& s);
template<> void scal<float>(size_t n, const float& a, float *x, int incx);
template<> void scal<double>(size_t n, const double& a, double *x, int incx);
template<> void scal< std::complex<float> >(size_t n, const std::complex<float>& a, std::complex<float> *x, int incx);
template<> void scal< std::complex<double> >(size_t n, const std::complex<double>& a, std::complex<double> *x, int incx);
template<> void scal<float>(size_t n, const float& a, std::complex<float> *x, int incx);
template<> void scal<double>(size_t n, const double& a, std::complex<double> *x, int incx);
template<> void swap<float>(size_t n, float *x, int incx, float *y, int incy);
template<> void swap<double>(size_t n, double *x, int incx, double *y, int incy);
template<> void swap< std::complex<float> >(size_t n, std::complex<float> *x, int incx, std::complex<float> *y, int incy);
template<> void swap< std::complex<double> >(size_t n, std::complex<double> *x, int incx, std::complex<double> *y, int incy);
template<> ptrdiff_t amax<float>(size_t n, const float *x, int incx);
template<> ptrdiff_t amax<double>(size_t n, const double *x, int incx);
template<> ptrdiff_t amax<float>(size_t n, const std::complex<float> *x, int incx);
template<> ptrdiff_t amax<double>(size_t n, const std::complex<double> *x, int incx);
// The accelerate framework doesn't appear to provide these
#ifndef USE_ACCELERATE
template<> ptrdiff_t amin<float>(size_t n, const float *x, int incx);
template<> ptrdiff_t amin<double>(size_t n, const double *x, int incx);
template<> ptrdiff_t amin<float>(size_t n, const std::complex<float> *x, int incx);
template<> ptrdiff_t amin<double>(size_t n, const std::complex<double> *x, int incx);
#endif
template<> void gemv<float>( size_t m, size_t n, const float* a, int lda, const float* x, int incx, float* y, int incy, const float& alpha, const float& beta, trans_e trans);
template<> void gemv<double>( size_t m, size_t n, const double* a, int lda, const double* x, int incx, double* y, int incy, const double& alpha, const double& beta, trans_e trans);
template<> void gemv< std::complex<float> >( size_t m, size_t n, const std::complex<float>* a, int lda, const std::complex<float>* x, int incx, std::complex<float>* y, int incy, const std::complex<float>& alpha, const std::complex<float>& beta, trans_e trans);
template<> void gemv< std::complex<double> >( size_t m, size_t n, const std::complex<double>* a, int lda, const std::complex<double>* x, int incx, std::complex<double>* y, int incy, const std::complex<double>& alpha, const std::complex<double>& beta, trans_e trans);
template<> void ger<float>(size_t m, size_t n, const float* x, int incx, const float* y, int incy, float* a, int lda, const float& alpha);
template<> void ger<double>(size_t m, size_t n, const double* x, int incx, const double* y, int incy, double* a, int lda, const double& alpha);
template<> void ger< std::complex<float> >(size_t m, size_t n, const std::complex<float>* x, int incx, const std::complex<float>* y, int incy, std::complex<float>* a, int lda, const std::complex<float>& alpha);
template<> void ger< std::complex<double> >(size_t m, size_t n, const std::complex<double>* x, int incx, const std::complex<double>* y, int incy, std::complex<double>* a, int lda, const std::complex<double>& alpha);
template<> void gerc<float>(size_t m, size_t n, const std::complex<float>* x, int incx, const std::complex<float>* y, int incy, std::complex<float>* a, int lda, const std::complex<float>& alpha);
template<> void gerc<double>(size_t m, size_t n, const std::complex<double>* x, int incx, const std::complex<double>* y, int incy, std::complex<double>* a, int lda, const std::complex<double>& alpha);
template<> void symv<float>(size_t n, const float* a, int lda, const float* x, int incx, float* y, int incy, const float& alpha, const float& beta, uplo_e upper);
template<> void symv<double>(size_t n, const double* a, int lda, const double* x, int incx, double* y, int incy, const double& alpha, const double& beta, uplo_e upper);
template<> void syr<float>(size_t n, const float *x, int incx, float *a, int lda, const float& alpha, uplo_e upper);
template<> void syr<double>(size_t n, const double *x, int incx, double *a, int lda, const double& alpha, uplo_e upper);
template<> void syr2<float>(size_t n, const float *x, int incx, const float *y, int incy, float *a, int lda, const float& alpha, uplo_e upper);
template<> void syr2<double>(size_t n, const double *x, int incx, const double *y, int incy, double *a, int lda, const double& alpha, uplo_e upper);
template<> void gemm<float>(size_t m, size_t n, size_t k, const float *a, int lda, const float *b, int ldb, float *c, int ldc, const float& alpha, const float& beta, conjtrans_e transa, conjtrans_e transb);
template<> void gemm<double>(size_t m, size_t n, size_t k, const double *a, int lda, const double *b, int ldb, double *c, int ldc, const double& alpha, const double& beta, conjtrans_e transa, conjtrans_e transb);
template<> void gemm< std::complex<float> >(size_t m, size_t n, size_t k, const std::complex<float> *a, int lda, const std::complex<float> *b, int ldb, std::complex<float> *c, int ldc, const std::complex<float>& alpha, const std::complex<float>& beta, conjtrans_e transa, conjtrans_e transb);
template<> void gemm< std::complex<double> >(size_t m, size_t n, size_t k, const std::complex<double> *a, int lda, const std::complex<double> *b, int ldb, std::complex<double> *c, int ldc, const std::complex<double>& alpha, const std::complex<double>& beta, conjtrans_e transa, conjtrans_e transb);
template<> void syrk<float>(size_t n, size_t k, const float *a, int lda, float *c, int ldc, const float& alpha, const float& beta, trans_e trans, uplo_e upper);
template<> void syrk<double>(size_t n, size_t k, const double *a, int lda, double *c, int ldc, const double& alpha, const double& beta, trans_e trans, uplo_e upper);
template<> void syrk< std::complex<float> >(size_t n, size_t k, const std::complex<float> *a, int lda, std::complex<float> *c, int ldc, const std::complex<float>& alpha, const std::complex<float>& beta, trans_e trans, uplo_e upper);
template<> void syrk< std::complex<double> >(size_t n, size_t k, const std::complex<double> *a, int lda, std::complex<double> *c, int ldc, const std::complex<double>& alpha, const std::complex<double>& beta, trans_e trans, uplo_e upper);
template<> void syr2k<float>(size_t n, size_t k, const float *a, int lda,const float *b, int ldb, float *c, int ldc, const float& alpha, const float& beta, trans_e trans, uplo_e upper);
template<> void syr2k<double>(size_t n, size_t k, const double *a, int lda,const double *b, int ldb, double *c, int ldc, const double& alpha, const double& beta, trans_e trans, uplo_e upper);
template<> void syr2k< std::complex<float> >(size_t n, size_t k, const std::complex<float> *a, int lda,const std::complex<float> *b, int ldb, std::complex<float> *c, int ldc, const std::complex<float>& alpha, const std::complex<float>& beta, trans_e trans, uplo_e upper);
template<> void syr2k< std::complex<double> >(size_t n, size_t k, const std::complex<double> *a, int lda,const std::complex<double> *b, int ldb, std::complex<double> *c, int ldc, const std::complex<double>& alpha, const std::complex<double>& beta, trans_e trans, uplo_e upper);
#endif


} // namespace blas
} // namespace warp

#endif
